Repositorio

Los diferentes análisis, preparación de datos y demás tareas requeridas en este trabajo fueron realizados utilizando R y Python. La información necesaria para replicar los resultados presentados en este reporte se encuentran en el siguiente repositorio de GitHub: https://github.com/psaldar/metodos-estadisticos-avanzados. Dentro de la información se incluyen los archivos que contienen tanto los datos antes y después de su procesamiento, además de los códigos de Python en Notebooks de Jupyter y el código de R en un documento como Markdown.

Introducción

El modelamiento del comportamiento de diferentes variables es un tema que ha sido estudiado en sectores energéticos, industriales, económicos y financieros. De allí se comienza a apreciar tanto la importancia que tienen los datos hoy en día al igual que las técnicas utilizadas para su modelamiento. La estadística es una disciplina que se preocupa por la recolección, organización, interpretación y análisis de datos, la cual, según su aplicación, puede traer un gran impacto en la industria al momento de la toma de decisiones. En particular, diferentes técnicas estadísticas han sido utilizadas en el sector financiero, las cuales permiten modelar comportamiento de los clientes, las acciones, y de otras variables.

Diferentes industrias dentro de su funcionamiento deben presentar ante la Superintendencia información relacionando los ingresos, los costos de venta (dinero invertido para la producción de un bien o servicio) y los gastos de venta (dinero invertido para la distribución de un bien o servicio) que presentaron anualmente. Además, se presume que del mercado colombiano es posible pensar que exista una relación entre las diferentes variables macroeconómicas (ej. PIB, TRM, Balance Fiscal, Indice de Desempleo, etc.) y estos montos de costos y gastos de las empresas. Debido a la cantidad de información con la que se cuenta (información de costos y gastos para diferentes empresas en colombia entre los años 2016 y 2018), se sabe que no se cuenta con información suficiente para la construcción de un modelo por empresa que permita ver la relación existente entre las variables macroeconómicas y las variables asociadas a costos y gastos de ventas. Por lo anterior, es posible considerar un conjunto de datos como la consolidación de la información de todas las empresas junto con la información macroeconomica para los años en estudio, así buscando construir un modelo general para realizar la modelación de estas variables reportadas ante la superintendencia para un conjunto de empresas cuya industria sea similar.

De esta manera, para el conjunto de datos mencionado anteriormente, nos enfocaremos en modelar la información de costos de venta a partir de las variables macroeconómicas disponibles. Se utilizarán modelos lineales para este trabajo: se usa un modelo lineal general y un modelo de efectos mixtos. Esta última estructura de modelos es bastante usada al momento de tener individuos que comparten la misma información pero tienen una salida diferente (en nuestro caso, todas las empresas comparten la misma información de las variables macroeconómicas pero tienen distintos costos de venta), por lo que utilizar este tipo de modelos resulta de gran interés ya que permite modelar tanto efectos a nivel de individuo como a nivel poblacional general.

Clasificación Industrial Seleccionada

El sector de la construcción es uno de los más relevantes en la economía colombiana. En nuestro contexto nacional el sector es considerado como uno de los más vitales para el desarrollo del país y representa uno de los más importantes rubros en materia de produccción interna componiendo cada año de 6 a 7 por ciento del producto interno bruto total y hasta un 7.1% del total de ocupados a nivel nacional. Dicho sector es caracterizado por sus fluctuaciones estacionarias fuertemente influenciadas por los planes de infraestructura de gran escala y los planes de gobierno. Respecto al último trimestre de 2019 tuvo un aumento de 3.4%, uno de los más altos a nivel de América latina.

Si bien las estimaciones globales para el año 2020 en Colombia para esta industria eran positivas, la coyuntura del COVID-19 ha de perturbar fuertemente el sector, afectando con alto impacto a los importadores de materiales y a la demanda frente a los retrasos en la ejecución de obras. Desde enero el sector de la producción de concreto ya estaba presentando caídas significativas de hasta un 8.3% con respecto al año pasado en el mismo mes, por lo que se esperan peores resultados al cierre del segundo trimestre del año presente. Muchas compañías planearon incrementos en sus precios, con la esperanza de generar mayores ingresos, pero dichos planes han de ser postergados bajo la actual coyuntura. La demanda, el driver más relevante en la industria, claramente se ve desplazado por efecto del impedimento de comercialización y la paralisis en el país bajo las medidas de cuarentena nacional. Un punto importante es que algunas compañías del sector podrán seguir con sus actividades producto de la inclusión de actividades de infraestructura como vital durante la crisis del COVID-19.

Frente a la incertidumbre que generan estos escenarios y teniendo en cuenta que la dinámica particular de cada compañía que aporta al crecimiento de la industria total evoluciona en el tiempo, se torna relevante construir análisis y modelos estadísticos robustos que respondan a las perturbaciones en el mercado y permitan obtener información valiosa para la toma de decisiones.

Consolidación de los datos (variables de salida y explicativas)

NOTA: toda esta parte inicial de preparación y consolidación de los datos se hizo en Python en un notebook de Jupyter. El notebook que realiza esta consolidación se encuentra en el repositorio, en la carpeta “PreparacionDatos”. Este script de R toma solo el dataset final que fue resultado de toda la consolidación y preprocesamiento.

Los datos de la variable de respuesta (costos y gastos de ventas de algunas empresas para los años 2016 a 2018) se extrajeron de: http://pie.supersociedades.gov.co.

Se conservaron solo las empresas cuya clasificación industrial uniforme internacional tenga la palabra “construcción” en ella y cuyos NIT estuvieran incluidos en el archivo de empresas de 2018 (Plenas e Individuales). Luego de aplicar este filtro, se eliminaron también todas aquellas empresas que en algún año entre 2015 y 2018 tuvieran algún valor vacío o nulo (NaN o cero) ya sea en sus costos de venta o en sus gastos de venta (se asumen que costos o gastos de venta de cero serían un error en la medición), dejando en total 33 empresas distintas en el dataset. A pesar de que solo se analizan en este reporte los años 2016 a 2018, fue necesario también descargar 2015 para poder calcular la variación porcentual anual para 2016. Para calcular las variables de costos y gastos de ventas como variación porcentual anual, s usa la siguiente fórmula:

\[ CostoDeVentasDif(t) = \dfrac{CostoDeVentas(t) – CostoDeVentas(t-1)}{|CostoDeVentas(t-1)|} \]

Donde t era el año. Así, se obtuvieron Costos_de_ventas_dif para los años 2016, 2017 y 2018. Se incluye el valor absoluto en el denominador en la fórmula para casos donde la variable sea negativa. Esto no ocurre con costos de ventas claramente, pero algunas variables macroeconómicas si tienen valores negativos como por ejemplo el Balance_Fiscal, por eso se generaliza la fórmula de variación porcentual para todas las variables con el uso del valor absoluto.

NOTA: A pesar de que en este trabajo los modelos que usamos son solo para el costo de ventas, aplicamos la misma transformación para el gasto de ventas para los análisis descriptivos.

Los valores para las variables macroeconómicas (variables explicativas) se extrajeron del Banco de la República y del DANE. Se descargaron archivos de Excel para cada una de las variables macro. Algunas de estas variables macro ya venían como variación porcentual entre un año y el otro. Sin embargo, otras no venían de esta manera y fue necesario transformarlas algunas de éstas. A continuación, detallamos la forma final en la que queda cada una de las variables (luego de nuestro preprocesamiento interno en Python).

Se hacen un par de aclaraciones:

Lectura de los datos

Aquí se lee el archivo de datos que ya fue preparado en Python, tal como se mencionó anteriormente.

all_data = read.csv("PreparacionDatos/Datos_completos.csv", encoding = "UTF-8")
NIT = all_data['NIT']

Transformación de los Datos

Para fines de intepretación de los modelos creados en el trabajo y para estabilizar un poco la varianza, consideraremos la transformación logarítmica tanto para los predictores como para la variable respuesta para poder tener un modelo de tipo log-log. Esto permitirá poder interpretar un coeficiente (beta) de la regresión como el incremento porcentual que se tendría de la variable dependiente al incrementar un 1% de la variable independiente relacionada con este coeficiente. Esta transformación del logaritmo también es de utilidad ya que se trabaja con porcentajes y en estos casos es usual que se presenten sesgos en la variable. De esta forma, la transformación logarítimica que se realiza es:

\[T(x) = log(1+x) \]

Por ejemplo, la variable de salida tendra la forma

\[ T(CostoDeVentasDif) = log(1+\dfrac{CostoDeVentas(t) – CostoDeVentas(t-1)}{|CostoDeVentas(t-1)|})\]

### Transformacion de las variables

### Transformacion a nivel logaritmico
all_data$Costo.de.ventas_dif_sinlog = all_data$Costo.de.ventas_dif
all_data$Costo.de.ventas_dif = log(1 + all_data$Costo.de.ventas_dif)
all_data$Gastos.de.ventas_dif_sinlog = all_data$Gastos.de.ventas_dif
all_data$Gastos.de.ventas_dif = log(1 + all_data$Gastos.de.ventas_dif)
all_data$PIB = log(1 + all_data$PIB)
all_data$TRM = log(1 + all_data$TRM)
all_data$Desempleo = log(1 + all_data$Desempleo)
all_data$Balance_CC = log(1 + all_data$Balance_CC)
all_data$Balance_Fiscal = log(1 + all_data$Balance_Fiscal)
all_data$Inflacion = log(1 + all_data$Inflacion)
all_data$Tasa_Intervencion = log(1 + all_data$Tasa_Intervencion)

Cabe mencionar que al momento de realizar las interpretaciones de los resultados del modelo en cuanto a su ajuste y su capacidad predictiva (en los plots que comparan valores predichos contra reales), se regresa la variable de salida a su representación de variaciones porcentuales anuales.

Análisis descriptivo

Damos una visualizacion inicial al conjunto de datos con el que se va a trabajar. Las variables de salida con sufijo _dif representan la variación porcentual ya con la transformada del logaritmo, las variables de salida sin sufijo son las variables originales en niveles y las variables de salida con sufijo _dif_sinlog representan las variaciones porcentuales anuales (sin aplicar aún la transformada del logaritmo). También se recuerda que ya se aplicó la transformación del logaritmo para todas la variables explicativas presentadas aquí y en las siguientes secciones.

Tener en cuenta que en la tabla siguiente, las primeras 33 observaciones son para el año 2016, las siguientes 33 para el año 2017 y las últimas 33 para el año 2018.

data = all_data[, -c(1,2)]

ref_ds=round( all_data[, -c(2)],4)
col_order_data <- c("NIT", "Costo.de.ventas", "Gastos.de.ventas", "Costo.de.ventas_dif", "Gastos.de.ventas_dif", "Costo.de.ventas_dif_sinlog", "Gastos.de.ventas_dif_sinlog", "TRM" ,"PIB", "Desempleo",  "Inflacion", "Tasa_Intervencion", "Balance_CC", "Balance_Fiscal")
ref_ds <- ref_ds[, col_order_data]

ref_ds$Costo.de.ventas <- currency(ref_ds$Costo.de.ventas, digits = 0L)
ref_ds$Gastos.de.ventas  <- currency(ref_ds$Gastos.de.ventas, digits = 0L)
#ref_ds$Costo.de.ventas_dif_sinlog<- percent(ref_ds$Costo.de.ventas_dif_sinlog, format = "f", digits = 4L)
#ref_ds$Gastos.de.ventas_dif_sinlog <- percent(ref_ds$Gastos.de.ventas_dif_sinlog   , format = "f", digits = 4L)
#ref_ds$TRM<- percent(ref_ds$TRM, format = "f", digits = 4L)
#ref_ds$PIB<- percent(ref_ds$PIB, format = "f", digits = 4L)
#ref_ds$Desempleo<- percent(ref_ds$Desempleo, format = "f", digits = 4L)
#ref_ds$Inflacion<- percent(ref_ds$Inflacion, format = "f", digits = 4L)
#ref_ds$Tasa_Intervencion<- percent(ref_ds$Tasa_Intervencion, format = "f", digits = 4L)
#ref_ds$Balance_CC<- percent(ref_ds$Balance_CC, format = "f", digits = 4L)
#ref_ds$Balance_Fiscal<- percent(ref_ds$Balance_Fiscal, format = "f", digits = 4L)


kable(ref_ds, digits = 4) %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T, font_size = 12)%>%
  scroll_box("100%", height = "250px")
NIT Costo.de.ventas Gastos.de.ventas Costo.de.ventas_dif Gastos.de.ventas_dif Costo.de.ventas_dif_sinlog Gastos.de.ventas_dif_sinlog TRM PIB Desempleo Inflacion Tasa_Intervencion Balance_CC Balance_Fiscal
800015615 $60,626,083 $13,978,834 -0.1040 0.1307 -0.0988 0.1397 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
800045720 $3,785,517 $60,608 0.0716 0.2824 0.0742 0.3263 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
800081030 $39,634,059 $5,645,885 0.2912 0.5275 0.3380 0.6947 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
800112440 $24,304,534 $8,977,594 -0.0547 -0.0369 -0.0532 -0.0362 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
800118660 $31,280,967 $6,702,215 -0.0480 -0.1945 -0.0469 -0.1768 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
800157469 $11,886,427 $1,050,755 0.1186 0.0248 0.1259 0.0251 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
800232356 $104,518,885 $9,540,000 0.0311 0.1862 0.0315 0.2047 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
800236890 $23,301,270 $404,719 0.7364 0.0814 1.0884 0.0848 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
801002644 $94,444,630 $27,857,644 0.2877 0.1753 0.3333 0.1915 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
805012368 $42,291,070 $7,314,503 -0.0169 0.1538 -0.0168 0.1662 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
806014553 $54,866,113 $17,048,833 -0.3100 -0.2376 -0.2665 -0.2115 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
830030574 $13,969,323 $870,261 -0.1243 -0.6014 -0.1169 -0.4520 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
830037495 $43,762,061 $16,325,119 0.1738 -0.0236 0.1898 -0.0233 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
830052054 $61,208,106 $1,865,920 -0.1948 -0.1148 -0.1770 -0.1085 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
860009694 $187,349,222 $7,091,362 -0.0913 -0.1789 -0.0872 -0.1638 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
860030360 $17,821,805 $4,121,173 -0.0603 -0.1769 -0.0585 -0.1621 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
860033653 $31,824,537 $8,364,174 0.0692 -0.0976 0.0717 -0.0930 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
860050956 $40,609,990 $5,978,160 -0.0529 -0.0342 -0.0515 -0.0336 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
860501682 $39,165,473 $5,786,553 -0.1198 -0.0091 -0.1129 -0.0091 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
890117431 $20,467,343 $1,920,789 0.2938 -0.0227 0.3415 -0.0224 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
890300012 $39,750,093 $3,060,095 0.1059 0.1395 0.1117 0.1497 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
890311366 $32,850,914 $13,043,626 -0.0226 0.1610 -0.0223 0.1746 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
890904459 $76,957,720 $12,679,101 0.1561 0.5171 0.1690 0.6771 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
890909034 $5,580,194 $1,594,136 -1.7471 0.2480 -0.8257 0.2814 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
890929951 $79,806,625 $1,523,941 0.1805 -0.1076 0.1979 -0.1020 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
900173460 $6,542,677 $762,613 -0.4232 -0.6234 -0.3450 -0.4639 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
900184722 $50,453,000 $5,204,078 0.0858 -0.1280 0.0896 -0.1201 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
900204182 $105,696,206 $9,676,876 0.0454 0.0708 0.0464 0.0734 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
900234565 $49,639,368 $4,707,006 -0.1160 -0.2243 -0.1095 -0.2009 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
900364670 $24,874,616 $553,453 -0.0081 0.0257 -0.0081 0.0260 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
900378893 $192,984,673 $8,658,834 -0.3404 0.0061 -0.2885 0.0061 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
900389088 $4,619,968 $102,497 -0.4814 0.2076 -0.3821 0.2308 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
900437650 $18,774,860 $4,641,379 -2.0388 0.3834 -0.8698 0.4673 0.1059 0.0709 0.0882 0.0559 0.0686 0.3013 0.2411
800015615 $59,080,156 $16,739,919 -0.0258 0.1803 -0.0255 0.1975 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
800045720 $3,977,633 $44,455 0.0495 -0.3099 0.0508 -0.2665 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
800081030 $35,278,433 $5,221,895 -0.1164 -0.0781 -0.1099 -0.0751 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
800112440 $23,453,024 $9,648,892 -0.0357 0.0721 -0.0350 0.0748 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
800118660 $33,462,570 $6,492,463 0.0674 -0.0318 0.0697 -0.0313 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
800157469 $10,452,870 $1,024,366 -0.1285 -0.0254 -0.1206 -0.0251 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
800232356 $104,940,884 $9,722,277 0.0040 0.0189 0.0040 0.0191 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
800236890 $8,321,068 $385,350 -1.0297 -0.0490 -0.6429 -0.0479 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
801002644 $113,667,103 $35,312,207 0.1853 0.2371 0.2035 0.2676 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
805012368 $46,330,434 $8,786,914 0.0912 0.1834 0.0955 0.2013 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
806014553 $57,657,912 $14,583,064 0.0496 -0.1562 0.0509 -0.1446 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
830030574 $13,700,631 $684,869 -0.0194 -0.2396 -0.0192 -0.2130 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
830037495 $28,235,186 $10,682,060 -0.4382 -0.4241 -0.3548 -0.3457 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
830052054 $41,704,463 $1,926,936 -0.3837 0.0322 -0.3186 0.0327 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
860009694 $175,304,403 $7,528,803 -0.0665 0.0599 -0.0643 0.0617 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
860030360 $17,199,886 $4,543,229 -0.0355 0.0975 -0.0349 0.1024 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
860033653 $33,298,022 $8,454,121 0.0453 0.0107 0.0463 0.0108 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
860050956 $33,396,041 $4,565,071 -0.1956 -0.2697 -0.1776 -0.2364 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
860501682 $45,097,539 $6,602,777 0.1410 0.1320 0.1515 0.1411 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
890117431 $17,500,906 $2,139,796 -0.1566 0.1080 -0.1449 0.1140 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
890300012 $35,494,740 $2,887,493 -0.1132 -0.0581 -0.1071 -0.0564 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
890311366 $27,816,455 $11,027,936 -0.1664 -0.1679 -0.1533 -0.1545 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
890904459 $82,314,130 $9,064,583 0.0673 -0.3356 0.0696 -0.2851 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
890909034 $2,163,668 $874,776 -0.9474 -0.6001 -0.6123 -0.4513 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
890929951 $63,594,391 $1,087,341 -0.2271 -0.3376 -0.2031 -0.2865 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
900173460 $4,558,220 $548,247 -0.3614 -0.3300 -0.3033 -0.2811 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
900184722 $47,025,390 $4,787,114 -0.0704 -0.0835 -0.0679 -0.0801 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
900204182 $86,095,148 $8,509,424 -0.2051 -0.1286 -0.1854 -0.1206 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
900234565 $54,212,946 $4,813,622 0.0881 0.0224 0.0921 0.0227 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
900364670 $19,842,088 $373,644 -0.2260 -0.3929 -0.2023 -0.3249 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
900378893 $174,006,870 $6,260,025 -0.1035 -0.3244 -0.0983 -0.2770 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
900389088 $8,713,342 $163,706 0.6345 0.4682 0.8860 0.5972 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
900437650 $9,165,656 $4,647,942 -0.7171 0.0014 -0.5118 0.0014 -0.0341 0.0636 0.0896 0.0401 0.0595 0.1390 -0.2597
800015615 $66,347,258 $18,727,983 0.1160 0.1122 0.1230 0.1188 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
800045720 $3,715,164 $34,368 -0.0683 -0.2574 -0.0660 -0.2269 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
800081030 $38,817,595 $4,277,083 0.0956 -0.1996 0.1003 -0.1809 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
800112440 $28,140,401 $10,497,079 0.1822 0.0843 0.1999 0.0879 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
800118660 $36,901,472 $5,464,218 0.0978 -0.1724 0.1028 -0.1584 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
800157469 $12,068,388 $1,094,795 0.1437 0.0665 0.1546 0.0688 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
800232356 $116,765,248 $10,138,362 0.1068 0.0419 0.1127 0.0428 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
800236890 $591,366 $2,759,408 -2.6441 1.9686 -0.9289 6.1608 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
801002644 $107,841,149 $35,411,964 -0.0526 0.0028 -0.0513 0.0028 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
805012368 $57,190,250 $9,895,818 0.2106 0.1188 0.2344 0.1262 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
806014553 $72,711,967 $17,780,360 0.2320 0.1982 0.2611 0.2192 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
830030574 $13,320,540 $892,387 -0.0281 0.2647 -0.0277 0.3030 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
830037495 $27,018,838 $9,736,407 -0.0440 -0.0927 -0.0431 -0.0885 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
830052054 $41,660,420 $2,114,010 -0.0011 0.0927 -0.0011 0.0971 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
860009694 $211,914,084 $8,088,444 0.1897 0.0717 0.2088 0.0743 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
860030360 $16,852,833 $4,493,496 -0.0204 -0.0110 -0.0202 -0.0109 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
860033653 $35,255,510 $8,737,644 0.0571 0.0330 0.0588 0.0335 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
860050956 $30,663,434 $4,533,365 -0.0854 -0.0070 -0.0818 -0.0069 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
860501682 $57,703,425 $7,796,907 0.2465 0.1662 0.2795 0.1809 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
890117431 $18,458,595 $2,098,534 0.0533 -0.0195 0.0547 -0.0193 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
890300012 $24,601,978 $2,365,611 -0.3666 -0.1994 -0.3069 -0.1807 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
890311366 $28,242,513 $8,103,166 0.0152 -0.3082 0.0153 -0.2652 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
890904459 $92,360,857 $10,030,521 0.1152 0.1013 0.1221 0.1066 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
890909034 $1,736,156 $183,814 -0.2201 -1.5600 -0.1976 -0.7899 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
890929951 $75,929,483 $1,180,434 0.1773 0.0821 0.1940 0.0856 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
900173460 $4,461,062 $686,750 -0.0215 0.2252 -0.0213 0.2526 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
900184722 $46,291,017 $5,169,423 -0.0157 0.0768 -0.0156 0.0799 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
900204182 $87,856,734 $8,628,978 0.0203 0.0140 0.0205 0.0140 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
900234565 $48,454,520 $6,727,499 -0.1123 0.3348 -0.1062 0.3976 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
900364670 $13,653,296 $208,080 -0.3738 -0.5854 -0.3119 -0.4431 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
900378893 $169,567,783 $8,326,660 -0.0258 0.2853 -0.0255 0.3301 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
900389088 $6,880,597 $166,062 -0.2361 0.0143 -0.2103 0.0144 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633
900437650 $8,351,676 $5,112,417 -0.0930 0.0952 -0.0888 0.0999 0.0018 0.0687 0.0924 0.0313 0.0426 -0.3202 -0.0633

En la siguiente tabla, vemos información descriptiva de las variables que se considerarán en la etapa de modelamiento del trabajo. Esto nos permite tener una idea inicial de las características del conjunto de datos. Se recuerda que en los modelos nos enfocaremos solamente en la variable de salida Costo_de_ventas_dif y en cómo modelarla usando las variables explicativas. Sin embargo, se incluyen en esta parte del análisis explicativo otras variables de costos y gastos de venta.

col_order_data <- c("Costo.de.ventas", "Gastos.de.ventas", "Costo.de.ventas_dif", "Gastos.de.ventas_dif", "Costo.de.ventas_dif_sinlog", "Gastos.de.ventas_dif_sinlog", "TRM" ,"PIB", "Desempleo",  "Inflacion", "Tasa_Intervencion", "Balance_CC", "Balance_Fiscal")

data_summ <- data[, col_order_data]

summary_df = stat.desc(data_summ)

summary_df$Costo.de.ventas <- currency(summary_df$Costo.de.ventas, digits = 0L)
summary_df$Gastos.de.ventas  <- currency(summary_df$Gastos.de.ventas, digits = 0L)
#summary_df$Costo.de.ventas_dif_sinlog<- percent(summary_df$Costo.de.ventas_dif_sinlog, format = "f", digits = 4L)
#summary_df$Gastos.de.ventas_dif_sinlog <- percent(summary_df$Gastos.de.ventas_dif_sinlog   , format = "f", digits = 4L)
#summary_df$TRM<- percent(summary_df$TRM, format = "f", digits = 4L)
##summary_df$PIB<- percent(summary_df$PIB, format = "f", digits = 4L)
#summary_df$Desempleo<- percent(summary_df$Desempleo, format = "f", digits = 4L)
#summary_df$Inflacion<- percent(summary_df$Inflacion, format = "f", digits = 4L)
#summary_df$Tasa_Intervencion<- percent(summary_df$Tasa_Intervencion, format = "f", digits = 4L)
#summary_df$Balance_CC<- percent(summary_df$Balance_CC, format = "f", digits = 4L)
#summary_df$Balance_Fiscal<- percent(summary_df$Balance_Fiscal, format = "f", digits = 4L)

summary_df <- summary_df[c('min', 'max', 'range', 'median', 'mean', 'std.dev'),]

kable(summary_df, digits =4) %>%
  kable_styling(bootstrap_options = "striped", full_width = F, 
                position = "center", fixed_thead = T)%>%
  scroll_box("100%", height = "300px")
Costo.de.ventas Gastos.de.ventas Costo.de.ventas_dif Gastos.de.ventas_dif Costo.de.ventas_dif_sinlog Gastos.de.ventas_dif_sinlog TRM PIB Desempleo Inflacion Tasa_Intervencion Balance_CC Balance_Fiscal
min $591,366 $34,368 -2.6441 -1.5600 -0.9289 -0.7899 -0.0341 0.0636 0.0882 0.0313 0.0426 -0.3202 -0.2597
max $211,914,084 $35,411,964 0.7364 1.9686 1.0884 6.1608 0.1059 0.0709 0.0924 0.0559 0.0686 0.3013 0.2411
range $211,322,718 $35,377,596 3.3805 3.5287 2.0173 6.9507 0.1400 0.0073 0.0042 0.0246 0.0260 0.6215 0.5008
median $35,278,433 $5,204,078 -0.0258 0.0107 -0.0255 0.0108 0.0018 0.0687 0.0896 0.0401 0.0595 0.1390 -0.0633
mean $48,030,668 $6,552,627 -0.1051 -0.0118 -0.0380 0.0627 0.0246 0.0677 0.0901 0.0424 0.0569 0.0400 -0.0273
std.dev $45,447,021 $6,576,354 0.4434 0.3393 0.2783 0.6624 0.0597 0.0031 0.0017 0.0102 0.0108 0.2646 0.2071

Además de la información descriptiva presentada en la tabla anterior, podemos ver para cada una de las variables y su distribución de forma visual con la ayuda de la creación de histogramas de frecuencia.

Vamos a obtener la información relacionada a las medidas de centralidad y dispersión del conjunto de datos. Notar que en la tabla descriptiva presentada anteriormente, contamos con la presencia del vector de medias y medianas que describen la centralidad del conjunto de datos.

Boxplot para variables absolutas

par(mfrow=c(1,2))
bxplot_costos = boxplot(Costo.de.ventas~Year, data = all_data)
bxplot_gastos = boxplot(Gastos.de.ventas~Year, data = all_data)

### Outliers para costos
nits_tab = t(all_data[(all_data$Costo.de.ventas %in% bxplot_costos$out),]$NIT)
rownames(nits_tab) <- c("NITs")

kable(nits_tab)%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  column_spec(1,bold=T)
NITs 860009694 900378893 860009694 900378893 860009694 900378893 
### Outliers para gastos
nits_tab = t(all_data[(all_data$Gastos.de.ventas %in% bxplot_gastos$out),]$NIT)
rownames(nits_tab) <- c("NITs")

kable(nits_tab)%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  column_spec(1,bold=T)
NITs 801002644 801002644 800015615 801002644

Boxplot para variables originales (en niveles)

par(mfrow=c(1,2))
bxplot_costos = boxplot(Costo.de.ventas~Year, data = all_data)
bxplot_gastos = boxplot(Gastos.de.ventas~Year, data = all_data)

### Outliers para costos
nits_tab = t(all_data[(all_data$Costo.de.ventas %in% bxplot_costos$out),]$NIT)
rownames(nits_tab) <- c("NITs")

kable(nits_tab)%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  column_spec(1,bold=T)
NITs 860009694 900378893 860009694 900378893 860009694 900378893 
### Outliers para gastos
nits_tab = t(all_data[(all_data$Gastos.de.ventas %in% bxplot_gastos$out),]$NIT)
rownames(nits_tab) <- c("NITs")

kable(nits_tab)%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  column_spec(1,bold=T)
NITs 801002644 801002644 800015615 801002644

Boxplot para variaciones porcentuales transformadas con logaritmo

par(mfrow=c(1,2))
bxplot_costos_dif = boxplot(Costo.de.ventas_dif~Year, data = all_data)
bxplot_gastos_dif = boxplot(Gastos.de.ventas_dif~Year, data = all_data)

### Outliers para costos dif
nits_tab = t(all_data[(all_data$Costo.de.ventas_dif %in% bxplot_costos_dif$out),]$NIT)
rownames(nits_tab) <- c("NITs")

kable(nits_tab)%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  column_spec(1,bold=T)
NITs 800236890 890909034 900389088 900437650 800236890 890909034 900389088 900437650 800236890 890300012 900364670 
### Outliers para gastos dif
nits_tab = t(all_data[(all_data$Gastos.de.ventas_dif %in% bxplot_gastos_dif$out),]$NIT)
rownames(nits_tab) <- c("NITs")

kable(nits_tab)%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  column_spec(1,bold=T)
NITs 830030574 900173460 800045720 800236890 890311366 890909034 900234565 900364670

Boxplot para las variaciones porcentuales sin hacer aún la transformación logarítmica

par(mfrow=c(1,2))
bxplot_costos_dif_sinlog = boxplot(Costo.de.ventas_dif_sinlog~Year, data = all_data)
bxplot_gastos_dif_sinlog = boxplot(Gastos.de.ventas_dif_sinlog~Year, data = all_data)

### Outliers para costos dif sin log
nits_tab = t(all_data[(all_data$Costo.de.ventas_dif_sinlog %in% bxplot_costos_dif_sinlog$out),]$NIT)
rownames(nits_tab) <- c("NITs")

kable(nits_tab)%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  column_spec(1,bold=T)
NITs 800236890 890909034 900437650 800236890 890909034 900389088 800236890 
### Outliers para gastos dif sin log
nits_tab = t(all_data[(all_data$Gastos.de.ventas_dif_sinlog %in% bxplot_gastos_dif_sinlog$out),]$NIT)
rownames(nits_tab) <- c("NITs")

kable(nits_tab)%>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  column_spec(1,bold=T)
NITs 800081030 890904459 900389088 800045720 800236890 890311366 890909034 900234565 900364670 900378893

Análisis:

Al observar la distribución del comportamiento de las variables de ventas y gastos originales (en niveles) para los 3 años en consideración, podemos ver que no hay mucha variación entre el comportamiento de éstas a lo largo del tiempo, vemos que en particular, tanto el monto de costo y gasto promedio (mediana) es prácticamente constante para los 3 años, aunque se nota un aumento en los costos en el año 2018 para el caso de algunas empresas. Ademas, vemos que para el caso de los Costos, las mismas 2 empresas (relacionadas a los NITs 860009694 y 900378893) fueron las que presentaron un mayor valor de costos en los 3 periodos, por lo que el boxplot los considera registros outliers. Mirando la distribución de los gastos, notamos que una empresa tuvo los gastos mas altos en los 3 periodos (relacionadas al NIT 801002644) pero adicional, en el último año, la empresa con el NIT 800015615 tuvo un aumento notable en sus gastos, por lo que es considerado también como un outlier por el boxplot.

Pasando a observar los boxplots de la distribución de las variables porcentuales (transformadas y sin transformar) de costos y gastos, vemos para ambos casos los registros outliers son aquellos que tienen un valor tanto de disminución como incremento porcentual muy altos (nos damos cuenta de esto al ver los boxplots de la variación porcentual sin aplicar aún la transformación logarítimica). Podemos ver que para el caso del 2016 en la variable costos, las empresas outliers fueron las correspondientes a los NITs 800236890 y 890909034, de estas, la que no presenta una diferencia negativa fue la primera (800236890), por lo que puede indicar un aumento considerable (raro) en el mercado en término de los costos de venta (practicamente el doble del año anterior tal como se puede contrastar con la variable en niveles de costo de diferencia relativa), mientras que la otra empresa presenta una disminución considerable en este mismo concepto. Para el caso del 2017, una empresa fue considerada como outlier: la empresa considerada como registro raro (900389088) tuvo un aumento en sus costos de un 88% (casi el doble de los costos del año anterior. Este punto corresponde al punto mostrado en el boxplot de la variable de diferencia de costo de ventas sin log). Para el último periodo, las empresas con NITs 800045720, 800236890, 890311366, 890909034, 900234565, 900364670 y 900378893 fueron consideradas como outliers, solo se encuentra una de las empresas fue mencionada en uno de los dos periodos anteriores (800236890), la cual vuelve a reportar una disminución considerable en sus costos (esta empresa corresponde al punto mostrado por el boxplot de diferencia de costos de venta para el año 2018, contrastando así la disminución que presentó en este año).

Ahora mirando la variable de diferencia de gastos, solamente 2 empresas tienen comportamientos muy altos, presentando un aumento de casi el doble de los gastos anuales, estas empresas corresponden a los NITs 800081030 y 890904459 para los años 2016 y 2017 respectivamente, en cambio para el 2018, se tuvo variabilidad en los gastos de ventas de 8 empresas que son considerados anormales tanto por sus aumentos como disminuciones, aquellas empresas son las de los NIT 900389088, 800045720, 800236890, 890311366, 890909034, 900234565, 900364670, 900378893 y 900364670. En aquellas empresas que tuvieron aumentos considerables, reportan casi hasta un cuarto del aumento de gastos en comparación al año anterior, mientras que aquellas que reportan una disminución, vemos que sus gastos se disminuyeron hasta el doble del periodo pasado.

Así, estos cambios en las industrias en sus costos y gastos de ventas, pueden estar directamente relacionados con algunos eventos que hayan hecho la gran variación en la materia prima requerida para sus productos (tipos de mezlcas de cemento por ejemplo), además de que el sector seleccionado depende bastante de los planes de infraestructura que se tengan, por lo que también es factible que en los periodos analizados hayan ocurrido cambios en estos planes dentro de la contratación de cada empresa, y sean estas las causas que nos lleven a ver diferencias tan grandes de periodo a periodo (tanto positivas como negativas).

Nota

Ahora, habiendo ya analizado las distribuciones de las variables, pasamos entonces a analizar la relación de la variable respuesta con los predictores. De aquí en adelante, usaremos siempre la variable de respuesta transformada (tanto para estos análisis de relación con predictores como para los gráficos de contorno y posteriormente para todos los modelos), ya que la variable con la transformación logarítmica es la que usaremos para las regresiones log-log. De esta forma, en todas las siguientes secciones del reporte cuando se mencione la variable “Costo de ventas” es la variable original en niveles, y la variable “Costo_de_ventas_dif” es la transformación logarítmica de la variación porcentual anual de dicha variable. Lo mismo aplica para Gastos de Ventas y Gastos_de_ventas_dif (aunque estas de gastos de ventas solo se usan en el análisis descriptivo).

Miramos ahora las matrices de Covarianza y Correlación para tener una intuición de la variabilidad de la información que consideramos.

Matriz de Covarianzas

cov_data = cov(data[,col_order_data])
kable(formatC(cov_data,format = "e", digits = 2)) %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  scroll_box("100%", height = "480px")
Costo.de.ventas Gastos.de.ventas Costo.de.ventas_dif Gastos.de.ventas_dif Costo.de.ventas_dif_sinlog Gastos.de.ventas_dif_sinlog TRM PIB Desempleo Inflacion Tasa_Intervencion Balance_CC Balance_Fiscal
Costo.de.ventas 2.07e+15 1.35e+14 4.53e+06 1.09e+06 2.34e+06 -1.56e+06 7.87e+04 4.68e+03 -2.42e+02 6.62e+03 1.76e+03 2.72e+03 2.93e+05
Gastos.de.ventas 1.35e+14 4.32e+13 5.84e+05 3.53e+05 3.58e+05 1.04e+05 3.17e+03 3.39e+02 1.35e+02 -3.27e+02 -8.17e+02 -2.30e+04 1.49e+04
Costo.de.ventas_dif 4.53e+06 5.84e+05 1.97e-01 -3.60e-02 1.10e-01 -1.48e-01 6.71e-05 4.07e-05 3.50e-05 -1.39e-04 -2.15e-04 -5.64e-03 1.02e-03
Gastos.de.ventas_dif 1.09e+06 3.53e+05 -3.60e-02 1.15e-01 3.87e-03 1.86e-01 1.81e-03 1.42e-04 2.74e-05 1.73e-05 -1.62e-04 -5.20e-03 7.46e-03
Costo.de.ventas_dif_sinlog 2.34e+06 3.58e+05 1.10e-01 3.87e-03 7.75e-02 -3.65e-02 1.12e-03 9.12e-05 2.04e-05 -3.99e-06 -1.22e-04 -3.78e-03 4.68e-03
Gastos.de.ventas_dif_sinlog -1.56e+06 1.04e+05 -1.48e-01 1.86e-01 -3.65e-02 4.39e-01 8.77e-04 1.94e-04 1.34e-04 -4.85e-04 -8.16e-04 -2.18e-02 6.23e-03
TRM 7.87e+04 3.17e+03 6.71e-05 1.81e-03 1.12e-03 8.77e-04 3.56e-03 1.61e-04 -6.00e-05 5.01e-04 3.80e-04 7.96e-03 1.22e-02
PIB 4.68e+03 3.39e+02 4.07e-05 1.42e-04 9.12e-05 1.94e-04 1.61e-04 9.46e-06 -5.77e-07 1.39e-05 4.09e-06 1.88e-05 5.97e-04
Desempleo -2.42e+02 1.35e+02 3.50e-05 2.74e-05 2.04e-05 1.34e-04 -6.00e-05 -5.77e-07 3.05e-06 -1.68e-05 -1.89e-05 -4.60e-04 -1.62e-04
Inflacion 6.62e+03 -3.27e+02 -1.39e-04 1.73e-05 -3.99e-06 -4.85e-04 5.01e-04 1.39e-05 -1.68e-05 1.05e-04 1.05e-04 2.46e-03 1.54e-03
Tasa_Intervencion 1.76e+03 -8.17e+02 -2.15e-04 -1.62e-04 -1.22e-04 -8.16e-04 3.80e-04 4.09e-06 -1.89e-05 1.05e-04 1.17e-04 2.85e-03 1.03e-03
Balance_CC 2.72e+03 -2.30e+04 -5.64e-03 -5.20e-03 -3.78e-03 -2.18e-02 7.96e-03 1.88e-05 -4.60e-04 2.46e-03 2.85e-03 7.00e-02 2.02e-02
Balance_Fiscal 2.93e+05 1.49e+04 1.02e-03 7.46e-03 4.68e-03 6.23e-03 1.22e-02 5.97e-04 -1.62e-04 1.54e-03 1.03e-03 2.02e-02 4.29e-02

Matriz de Correlación

cor_data = cor(data[,col_order_data])
kable(round(cor_data,2)) %>%
  kable_styling(bootstrap_options = "striped", full_width = F, position = "center", fixed_thead = T)%>%
  scroll_box("100%", height = "480px")
Costo.de.ventas Gastos.de.ventas Costo.de.ventas_dif Gastos.de.ventas_dif Costo.de.ventas_dif_sinlog Gastos.de.ventas_dif_sinlog TRM PIB Desempleo Inflacion Tasa_Intervencion Balance_CC Balance_Fiscal
Costo.de.ventas 1.00 0.45 0.22 0.07 0.18 -0.05 0.03 0.03 0.00 0.01 0.00 0.00 0.03
Gastos.de.ventas 0.45 1.00 0.20 0.16 0.20 0.02 0.01 0.02 0.01 0.00 -0.01 -0.01 0.01
Costo.de.ventas_dif 0.22 0.20 1.00 -0.24 0.89 -0.50 0.00 0.03 0.05 -0.03 -0.04 -0.05 0.01
Gastos.de.ventas_dif 0.07 0.16 -0.24 1.00 0.04 0.83 0.09 0.14 0.05 0.00 -0.04 -0.06 0.11
Costo.de.ventas_dif_sinlog 0.18 0.20 0.89 0.04 1.00 -0.20 0.07 0.11 0.04 0.00 -0.04 -0.05 0.08
Gastos.de.ventas_dif_sinlog -0.05 0.02 -0.50 0.83 -0.20 1.00 0.02 0.10 0.12 -0.07 -0.11 -0.12 0.05
TRM 0.03 0.01 0.00 0.09 0.07 0.02 1.00 0.87 -0.58 0.82 0.59 0.50 0.99
PIB 0.03 0.02 0.03 0.14 0.11 0.10 0.87 1.00 -0.11 0.44 0.12 0.02 0.94
Desempleo 0.00 0.01 0.05 0.05 0.04 0.12 -0.58 -0.11 1.00 -0.94 -1.00 -1.00 -0.45
Inflacion 0.01 0.00 -0.03 0.00 0.00 -0.07 0.82 0.44 -0.94 1.00 0.95 0.91 0.73
Tasa_Intervencion 0.00 -0.01 -0.04 -0.04 -0.04 -0.11 0.59 0.12 -1.00 0.95 1.00 1.00 0.46
Balance_CC 0.00 -0.01 -0.05 -0.06 -0.05 -0.12 0.50 0.02 -1.00 0.91 1.00 1.00 0.37
Balance_Fiscal 0.03 0.01 0.01 0.11 0.08 0.05 0.99 0.94 -0.45 0.73 0.46 0.37 1.00 
corrplot(cor_data, method="circle")

Dado que en este trabajo se pretende realizar la modelación de alguna de las variables de costo de ventas (también considerando la transformación de las variables por periodo), la matriz de correlación es un buen indicador para observar relaciones lineales existentes en las variables. En particular, vemos correlaciones relativamente fuertes entre algunas de las versiones de las variables de costos y gastos. Para este caso, es natural encontrar una correlación alta en varias de estas variables, pues es usual que los costos de venta y los gastos de venta se relacionen entre sí, ya que ambas tienen que ver con dinero que las empresas invierten ya sea para producir un producto (costos de venta) o para distribuirlo (gastos de ventas). Además, es normal también ver una alta correlación entre una variable y su propia transformación.

Adicionalmente, se nota también que las variables macroeconómicas suelen tener altas correlaciones entre sí, esto era algo de esperar ya que todas están relacionadas con la economía del país y es de esperar que sí la economía esté mejorando ciertas variables incrementen y otras disminuyan, y viceversa.

Por otro lado, para la variable en la que nos enfocamos que es Costo_de_ventas_dif, no se ven correlaciones fuertes en relación a las variables macro-económicas, lo que nos indica que no existe una clara relación lineal entre esta variable y las demás. De esta matriz de correlación es claro que hay las dependencias lineales entre las variables explicativas y las variables de salida son muy débiles, lo que complica que haya un buen desempeño de los modelos. Sin embargo, esto era de esperar debido a que se tienen solo 3 valores únicos de las variables macroeconómicas mientras que las variables de respuesta tienen una gran variabilidad, por lo que muchas veces cambian los valores de la variable de respuesta y no cambia el valor de las variables explicativas, lo que conlleva a que sea baja su correlación.

Visualización de la distribución de los datos

Realizamos ahora una visualización (gráficos de dispersión) para todos los pares de variables que usaremos en los modelos, es decir, para la variable de respuesta Costo_de_ventas_dif y para las variables macroeconómicas (todas ya con la transformación logarítmica).

vars = c("Costo.de.ventas_dif","TRM","PIB","Inflacion","Desempleo","Tasa_Intervencion", "Balance_CC", "Balance_Fiscal" )
pairs(data[,vars],labels = vars,font.labels = 2,diag.panel = panel.hist, lower.panel = panel.cor)

En el gráfico anterior, podemos ver la información de la distribución de las variables que se utilizan en los modelos, vemos entonces que existen correlaciones fuertes entre las variables macroeconómicas (tal como nos lo habia indicado la matriz de correlación), al mismo tiempo que se ven correlaciones muy débiles entre la variable de respuesta y los predictores. A partir de estos resultados, es claro que aquella variable que tiene una mayor variabilidad es Costo_de_ventas_dif, pues observamos que para las variables macroeconómicas sólo contamos con la información de 3 años, y los valores que toman cada año son los mismos todos la individuos (empresas) que estamos considerando en este trabajo, llevando así a que tengamos solo 3 valores únicos para cada variable macroeconómica, mientras que se tienen 99 valores únicos para la variable de respuesta.

Visualización de los datos según curvas de normalidad

### Funcion para encontrar los contornos
c_alpha = function(alpha, sigma, p){
  res = (2*pi)^(-p/2)*(det(sigma))^(-1/2)*exp(-1/2*qchisq(1-alpha, df = p))
  
  return(res)
}

grafica = function(data, name1, name2){

  data_aux = data[c(name1,name2)]
  names(data_aux) = c("y1","y2")
  
  cov_data = cov(data_aux)
  mean_data = mean=colMeans(data_aux)
  
  min_value1 = min(data[name1])
  max_value1 = max(data[name1])
  
  min_value2 = min(data[name2])
  max_value2 = max(data[name2])
  
  n = 100
  
  y1 = seq(min_value1, max_value1, length.out = n)
  y2 = seq(min_value2, max_value2, length.out = n)
  
  grid = expand.grid(y1,y2)
  grid$Z<-apply(grid,1,dmvnorm,mean = mean_data,sigma=cov_data)
  Z<-matrix(grid$Z,nrow=n,ncol=n)
  
  contornos = sapply(c(0.01, 0.05, 0.1), c_alpha, sigma = cov_data, 2)
  
  contour(y1,y2,Z,levels=contornos,labels=c("99%","95%","90%"),
          las=1)
  points(data_aux$y1,data_aux$y2)
  grid()  
  title(main = "Contornos de distribucion normal", xlab = name1, ylab = name2)
}
p1 = grafica(all_data, "Costo.de.ventas_dif", "PIB")

p2 = grafica(all_data, "Costo.de.ventas_dif", "Desempleo")

p3 = grafica(all_data, "Costo.de.ventas" , "Costo.de.ventas_dif")

p4 = grafica(all_data, "Costo.de.ventas" , "Costo.de.ventas_dif_sinlog")

Las gráficas anteriores nos permiten ver el comportamiento de la distribución de dos variables, además de tener como referencia los contornos de una distribución normal bivariada. En estos casos nos enfocamos en la variable de Costos de ventas y sus distintas transformaciones. De estos gráficos, podemos concluir que aquellos pares de variables graficados, no son candidatas para afirmar que pueden seguir una distribución normal, pues claramente tenemos muchos registros dentro del conjunto de datos cuyos valores caen por fuera de las regiones de normalidad (más de los esperados en una distribución normal). Se esperaría que solo aproximadamente el 1% de los valores caiga por fuera de la región del 99% y está cayendo por fuera un porcentaje mucho mayor al 1%.

De igual manera, al considerar dos variables macroeconómicas (el PIB y el Desempleo), justo por la característica de las mismas en el conjunto de datos (tener sólo valores para 3 años), no nos dan información acerca de la dispersión de la variable, por el contrario, tiene un efecto en las gráficas de contorno de enfocar los datos en 3 lineas (una para cada año). La variabilidad observada en los contornos que contienen variables macroeconómicas depende casi completamente de la variable de respuesta, ya que esta variable sí tiene una alta variabilidad. De igual forma, en estos gráficos de contorno que comparan un par de variables explicativas con la variable de respuesta, vemos que encontramos muchos puntos de los indivuduos que caen por fuera de la región de normalidad, por lo tanto tampoco podemos concluir normlidad bivariada entre estas variables.

Creación de modelos

Recordar que, tal como se mencionó al inicio, se realizó una transformación logarítmica tanto a las variables macroeconómicas como a la variable respuesta (variación porcentual anual de los costos de venta), por lo que al momento de la interpretación de los modelos, estaremos trabajando con regresiones del estilo log-log.

data_aux = select(all_data,select = -starts_with('NIT'))
data_aux = select(data_aux,select = -starts_with('Year'))

Realizamos la partición del conjunto de datos en dos conjuntos diferentes, el primer conjunto corresponde a los datos utilizados en la creación o entrenamiento del modelo y será usado para evaluar el ajuste, mientras que el segundo corresponde a aquellos datos que el modelo no ha visto, por lo que servirá para probar el desempeño del modelo en su capacidad predictiva.

La selección de ambos conjuntos se basa en la información de las empresas, pues se desea que el valor promedio de las variaciones porcentuales anuales del costo de ventas de las empresas que se considera en el conjunto de entrenamiento, sea similar al valor promedio de las variaciones porcentuales anuales del costo de ventas de las empresas en el conjunto de prueba. Partiendo de ese objetivo para la partición, se consideran 100000 combinaciones diferentes de las empresas (combinaciones donde el conjunto de entrenamiento tenga 22 empresas y el conjunto de prueba tenga 11 empresas), y se selecciona la combinación de empresas que tenga la menor diferencia al cuadrado entre los promedios de las variaciones porcentuales anuales del costo de ventas de las empresas de entrenamiento y de prueba. Esta tarea de la selección de las empresas fue realizada en un Notebook de python que se puede encontrar en el repositorio en la carpeta “PreparacionDatos”.

empresas_train = c(800015615, 800045720,
                   800081030, 800112440,
                   800118660, 800157469,
                   800232356, 800236890,
                   801002644, 805012368,
                   806014553, 830030574,
                   830037495, 860009694,
                   860033653, 860050956,
                   890909034, 900173460,
                   900184722, 900204182,
                   900364670, 900378893)

empresas_test = c(830052054, 860030360,
                   860501682, 890117431,
                   890300012, 890311366,
                   890904459, 890929951,
                   900234565, 900389088,
                   900437650)
train = all_data$NIT%in%empresas_train
test = all_data$NIT%in%empresas_test

NIT_train = NIT[train,]
NIT_test = NIT[test,]

data_train = data_aux[train,]
rownames(data_train) <- NULL
data_train_z = data_train


data_train_z = as.data.frame(data_train_z)

data_test = data_aux[test,]
rownames(data_test) <- NULL
data_test_z = data_test

data_train_z$NIT = NIT_train
data_test_z$NIT = NIT_test

Se crean funciones para calcular el R-squared y R-Squared Adjusted, las cuales serán usadas más adelante para evaluar el ajuste y la capacidad predictiva de los modelos.

r2_score <- function(x, y) summary(lm(y~x))$r.squared
adj_r2_score <- function(x, y) summary(lm(y~x))$adj.r.squared

Modelos para los Costos de venta

En esta sección, se usan un par de modelos (modelo lineal general y modelo lineal con efectos mixtos) para modelar la variable de Costos_de_venta_dif como una regresión de las variables macroeconómicas (todas ya transformadas, para que sea una regresión log-log). De igual forma, se aclara que cuando se creen los plots para medir ajuste y capacidad predictiva (los que comparan los valores predichos contra los reales), se aplica una transformación para poder retornar la variable de salida para que represente las variaciones porcentuales anuales (y no su logaritmo). Para esto, lo que se hace es aplicar la función exponencial a la variable de salida, y a ese valor se le resta 1, logrando así recuperar la escala de la variable para que represente una variación porcentual anual en los gráficos.

Modelo Lineal General

Primero, se evalúa el rendimiento de un modelo lineal general para realizar las predicciones. Se utiliza una selección stepwise para decidir que variables se incluirían en este modelo lineal.

step.model <- ols_step_both_p(lm('Costo.de.ventas_dif~PIB+TRM+Desempleo+Inflacion+Tasa_Intervencion+Balance_CC+Balance_Fiscal', data = data_train_z), pent = 0.1, prem = 0.1, details = TRUE)
## Stepwise Selection Method   
## ---------------------------
## 
## Candidate Terms: 
## 
## 1. PIB 
## 2. TRM 
## 3. Desempleo 
## 4. Inflacion 
## 5. Tasa_Intervencion 
## 6. Balance_CC 
## 7. Balance_Fiscal 
## 
## We are selecting variables based on p value...
## 
## 
## Stepwise Selection: Step 1 
## 
## - PIB added 
## 
##                           Model Summary                           
## -----------------------------------------------------------------
## R                        0.086       RMSE                  0.457 
## R-Squared                0.007       Coef. Var          -429.959 
## Adj. R-Squared          -0.008       MSE                   0.209 
## Pred R-Squared          -0.045       MAE                   0.240 
## -----------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.100         1          0.100    0.478    0.4919 
## Residual       13.367        64          0.209                    
## Total          13.467        65                                   
## ------------------------------------------------------------------
## 
##                                    Parameter Estimates                                    
## -----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig       lower     upper 
## -----------------------------------------------------------------------------------------
## (Intercept)    -0.967         1.246                 -0.776    0.441     -3.456     1.522 
##         PIB    12.709        18.384        0.086     0.691    0.492    -24.017    49.436 
## -----------------------------------------------------------------------------------------
## Warning in if (pvals[minp] <= pent) {: la condición tiene longitud > 1 y sólo el
## primer elemento será usado
## 
## No more variables to be added/removed.
## 
## 
## Final Model Output 
## ------------------
## 
##                           Model Summary                           
## -----------------------------------------------------------------
## R                        0.086       RMSE                  0.457 
## R-Squared                0.007       Coef. Var          -429.959 
## Adj. R-Squared          -0.008       MSE                   0.209 
## Pred R-Squared          -0.045       MAE                   0.240 
## -----------------------------------------------------------------
##  RMSE: Root Mean Square Error 
##  MSE: Mean Square Error 
##  MAE: Mean Absolute Error 
## 
##                               ANOVA                                
## ------------------------------------------------------------------
##                Sum of                                             
##               Squares        DF    Mean Square      F        Sig. 
## ------------------------------------------------------------------
## Regression      0.100         1          0.100    0.478    0.4919 
## Residual       13.367        64          0.209                    
## Total          13.467        65                                   
## ------------------------------------------------------------------
## 
##                                    Parameter Estimates                                    
## -----------------------------------------------------------------------------------------
##       model      Beta    Std. Error    Std. Beta      t        Sig       lower     upper 
## -----------------------------------------------------------------------------------------
## (Intercept)    -0.967         1.246                 -0.776    0.441     -3.456     1.522 
##         PIB    12.709        18.384        0.086     0.691    0.492    -24.017    49.436 
## -----------------------------------------------------------------------------------------

De los resultados anteriores, vemos que el modelo lineal general con mejor desempeño es aquel que solo considera la variable relacionada con el PIB además de un intercepto. Por lo tanto, procedemos a mirar el desempeño de un modelo lineal general entrenado considerando únicamente la variable del PIB y un intercepto.

mod_lin = lm('Costo.de.ventas_dif~PIB', data = data_train_z)
summary(mod_lin)
## 
## Call:
## lm(formula = "Costo.de.ventas_dif~PIB", data = data_train_z)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.55042 -0.01670  0.08328  0.19822  0.80263 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.9668     1.2460  -0.776    0.441
## PIB          12.7095    18.3843   0.691    0.492
## 
## Residual standard error: 0.457 on 64 degrees of freedom
## Multiple R-squared:  0.007412,   Adjusted R-squared:  -0.008097 
## F-statistic: 0.4779 on 1 and 64 DF,  p-value: 0.4919

En este caso el modelo está siendo exclusivamente dependiente de la dinámica del PIB, y se tiene que un aumento de un 1% en la variación porcentual del PIB anual conllevaría a un aumento de 12.71% sobre la variación porcentual de los costos de venta, lo cual es razonable, en particular porque, al percibir mayor crecimiento en la economía global, la industria naturalmente también se ve afectada al haber mayor producción y un mayor retorno de dinero; no obstante, las dinámicas de las compañías que componen la industria muchas veces se ven influenciada por el factor macro más global de manera diferentes, ello por su tamaño o tendencias, así que probablemente sea necesario modelar dichos efectos de una forma más robusta. Esta dinámica explicada también tiene sentido para el sector particular analizado (el de la construcción): un aumento de la economía global (PIB) está relacionado con un crecimiento global de las empresas y por ende también de forma general de las empresas de la construcción. Al ser más grandes las empresas es normal que tengan mayores ingresos (por ventas de materiales contrucción o por proyectos de construcción) y eso conlleva también a que tengan mayores costos.

NOTA: se hace una aclaración aquí para evitar confusiones. Cuando se dice que “aumento de un 1% en la variación porcentual del PIB anual conllevaría a un aumento de 12.71%”, el aumento del 1% hace referencia a aumentar 1% de la variación porcentual, no sumarle 1% a esta. Esto quiere decir entonces que, si por ejemplo se tiene una variación porcentual anual del 5%, el aumento del 1% mencionado haría referencia en aumentar 5% en un 1% de si mismo (5% + 5%/100). Esto no se debe confundir con aumentar 5% y sumarle un 1% y llegar a 6%, ya que este sería un caso distinto y no es el caso analizado en la regresión log-log.

Desempeño del Modelo Entrenado

Miremos el desempeño del modelo en el conjunto de entrenamiento (capacidad de ajuste):

preds = predict(mod_lin)

r2_model<-r2_score(preds, data_train_z$Costo.de.ventas_dif)
adj_r2_model<-adj_r2_score(preds, data_train_z$Costo.de.ventas_dif)
mse_model<-mean((data_train_z$Costo.de.ventas_dif - preds)^2)

#sub_tit = paste("R2", format(r2_model, digits=2, nsmall=2), 
#                "; R2_adj", format(adj_r2_model, digits=2, nsmall=2), 
#                sep = " ", collapse = NULL)

plot(lapply(preds, function(x) exp(x) - 1),  
     lapply(data_train_z$Costo.de.ventas_dif, function(x) exp(x) - 1), 
     xlab = 'Valor predicho', ylab = 'Valor real', 
     main='Modelo Lineal General (variación porcentual anual)')#, sub = sub_tit)
abline(a=0, b=1, lwd = 2, col = 'red')

Medidas Resultados
r2 train 0.0074123
r2 adj train -0.0080969
mse train 0.2025358

Como puede apreciarse en la gráfica anterior el comportamiento de las estimaciones es relativamente estático y no logra distribuir los puntos de información de forma dinámica. En particular, es claro como se obtuvo una segmentación en tres partes, lo cual es razonable y va en línea con la idea de crecimiento, estancamiento y decaida de los costos de venta de las compañías de la industria, es decir, si bien los valores predichos se encuentran desvíados de los reales, el sólo PIB y la propia variación de los costos permiten identificar algunas tendencias en la industria. Este comportamiento de solo predecir 3 escenarios diferentes tiene mucho sentido ya que se sabe que solo existen 3 valores diferentes para las variables explicativas en todo el conjunto de datos.

Adicionalmente, al ver las métricas obtenidas por este modelo en el conjunto de entrenamiento comprobamos también que el ajuste realizado es muy pobre: tuvo un R2 de 0.0074 (muy cercano a 0), indicando un pobre ajuste. Además su R2 ajustado fue también muy bajo (-0.0081). Su MSE fue de 0.2025, el cuál puede ser mejorado notablemente (cómo se verá más adelante con los modelos de efectos mixtos).

Observando ya las cifras predichas, vemos que las predicciones en el conjunto de ajuste oscilan desde -0.15 hasta -0.06, lo que muestra que están generando predicciones razonables (cercanas a la media al no tener mucha más capacidad de mejora). Son razonables estas predicciones en cuenta la media de (-0.03) y desviación estandar (0.27) de la variable de respuesta (sin transformación logarítmica), las cuáles fueron presentadas en el análisis descriptivo. Estas predicciones se puede decir que son acordes con el sentido común (predicen una disminución porcentual de costos entre 6% y 15%, el cuál es un valor factible).

Ahora, evaluemos el rendimiento del modelo en el conjunto de prueba (capacidad predictiva):

preds_test = predict(mod_lin, newdata = data_test_z)

r2_model<-r2_score(preds_test, data_test_z$Costo.de.ventas_dif)
adj_r2_model<-adj_r2_score(preds_test, data_test_z$Costo.de.ventas_dif)
mse_model<-mean((data_test_z$Costo.de.ventas_dif - preds_test)^2)

sub_tit = paste("R2", format(r2_model, digits=2, nsmall=2), 
                "; R2_adj", format(adj_r2_model, digits=2, nsmall=2), 
                sep = " ", collapse = NULL)

plot(lapply(preds_test, function(x) exp(x) - 1), 
     lapply(data_test_z$Costo.de.ventas_dif, function(x) exp(x) - 1), 
     xlab = 'Valor predicho', ylab = 'Valor Real', main='Modelo Lineal General (variación porcentual anual)')#, sub=sub_tit)
abline(a=0, b=1, lwd = 2, col = 'red')
Medidas Resultados
r2 test 0.0083219
r2 adj test -0.0236677
mse test 0.1802350

La aplicación sobre el conjunto de prueba mostrada genera resultados muy similares a aquellos sobre el de entrenamiento. La capacidad predictiva es también considerablemente baja (al igual que el ajuste), lo cual es esperado, debido a la separación natural de las tendencias de las empresas y a la falta de la modelación del efecto del indicador sobre la medida de estudio; no obstante, el modelo entrega valores que se mantienen en un rango conservador, es decir, no genera estimaciones imposibles, sino que hace predicciones que efectivamente pudieran alcanzarse bajo reducciones en las cuentas de costos de ventas.

Adicionalmente, al ver las métricas obtenidas por este modelo en el conjunto de prueba comprobamos también que su capacidad predictiva es muy pobre: el MSE es de 0.1802, mientras que el R2 y el R2 ajustado son muy bajos (0.0083 y -0.0237 respectivamente).

Observando ya las cifras predichas, vemos que las predicciones en el conjunto de prueba oscilan desde -0.15 hasta -0.06, lo que muestra que están generando predicciones razonables (cercanas a la media al no tener mucha más capacidad de mejora). Estas predicciones se puede decir que son acordes con el sentido común (predicen una disminución porcentual de costos entre 6% y 15%, el cuál es un valor factible). Es importante notar que las cifras predichas son exactamente iguales a las predichas en el conjunto de entrenamiento, ya que las variables explicativas que el modelo ve en train y en test son las mismas (solo los 3 valores únicos del PIB, los cuáles ve de forma repetida en varias ocasiones). Sin embargo, obviamente las métricas obtenidas son distintas ya que las cifras reales de las observaciones de los conjuntos de entrenamiento y prueba sí son distintas.

Los resultados obtenidos tanto en la capacidad de ajuste del modelo lineal como en su capacidad predictiva son los esperados: al tener variables explicativas exactamente iguales para todas las empresas del mismo año, es de esperar que no logré obtener buenos resultados ya que no tiene cómo modelar de forma particular a cada empresa (no se le incluyen variables indicadoras). Sin embargo, esperamos resolver más adelante este problema al introductir un modelo de efectos mixtos para la estimación, el cuál logrará detectar los factores particulares de cada empresa y separar mejor el componente poblacional del particular.

Gráficos de Evaluación de modelos

Se observan algunos de los gráficos de evaluación de modelos para mirar el rendimiento del modelo lineal.

par(mfrow=c(3,2))
mod_plot = plot(mod_lin, which = c(1:6))

De los gráficos anteriores, en particular analizando el valor de la distancia de Cook para las observaciones, notamos que las observaciones 17, 30, 52 son candidatas a ser registros atípicos en el conjunto de datos. Notar que los outliers corresponden a la empresa con NIT 890909034 para el 2016, y para la empresa 800236890 en los periodos del 2017 y 2018, información que coincide con algunos de los outliers obtenidos en el boxplot. En vista de su comportamiento raro en magnitud según el boxplot, no es de extrañarse que estos registros hayan sido aquellos con la distancia de Cook mas grande.

Modelos Lineales Generalizados (GLM)

Dentro de los modelos que fueron expuestos a lo largo del curso, se encuentra la familia de los modelos lineales generalizados, los cual nos permiten la creación de modelos teniendo en cuenta diferentes supuestos que se hacen sobre la variable respuesta que estamos considerando. En la elaboración de este trabajo, se está realiza un modelo de predicción sobre transformaciones logarítmica en la diferencia relativa de los costos de ventas que tienen las empresas, por lo que el dominio de la variable son los números reales (ya que puede tomar valores continuos tanto positivos como negativos, los cuales son considerados en las diferentes transformaciones que se plantean). Justo por la caracteristica de la variable respuesta que tenemos, no se considera pertinente realizar evaluacion de los modelos, esta conclusión está apoyada del análisis que se realiza de los diferentes modelos generalizados que podriamos considerar:

Por el análisis realizado, no se considera pertinente utilizar alguno de los otros modelos lineales generalizados para el modelamiento de la variable transformada de diferencia relativa de costos de venta.

Modelo de Efectos Mixtos

Debido a que se tenían muchas observaciones con exactamente las mismas variables explicativas y donde solo cambiaba el individuo (en este caso, la empresa), este problema era un problema ideal para que fuera abordado por modelos de efectos mixtos. Por dicha razón, decidimos usar un modelo de efectos mixtos que utilice las variables explicativas relacionadas con el PIB y el desempleo para modelar tanto los efectos generales de la población como los efectos aleatorios correspondientes a cada empresa en particular.

A continuación, se presenta la estimación del modelo de efectos mixtos mencionado:

mod_me = lmer('Costo.de.ventas_dif~Desempleo+PIB+(0+Desempleo|NIT)+(PIB|NIT)', data = data_train_z)
summary(mod_me)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Costo.de.ventas_dif ~ Desempleo + PIB + (0 + Desempleo | NIT) +  
##     (PIB | NIT)
##    Data: data_train_z
## 
## REML criterion at convergence: 68.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.2577 -0.0832  0.1859  0.3549  2.5374 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  NIT      Desempleo   3.576169 1.89108       
##  NIT.1    (Intercept) 0.004193 0.06475       
##           PIB         0.566302 0.75253  -1.00
##  Residual             0.182955 0.42773       
## Number of obs: 66, groups:  NIT, 22
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  -0.8466     3.0990  -0.273
## Desempleo    -1.2765    30.4740  -0.042
## PIB          12.6315    17.3072   0.730
## 
## Correlation of Fixed Effects:
##           (Intr) Desmpl
## Desempleo -0.926       
## PIB       -0.473  0.107

Para este modelo de efectos mixtos, observando primero la componente de los efectos fijos, vemos que un incremento en un 1% de la variación del PIB anual del pais, reflejaría un aumento en un 12.63% en la variación porcentual de los costos de venta anuales de las empresas. Por otro lado, al ver un aumento de un 1% en la tasa de desempleo, se disminuiría la variación porcentual de los costos de venta en un 1.2765%. Dichos resultados son completamente razonables en la realidad por el argumento previamente expuesto sobre el efecto del crecimiento del PIB en la variación de los costos: al presentarse mayor crecimiento económica, es usual que las empresas tengan en general mayores ingresos y esto suele conllevar a más costos. Por otro lado, una relación negativa entre el desempleo y la diferencia de costos de venta claramente tiene sentido ya que al presentarse mayor desempleo se tendrá una mano de obra más barata lo que reduce los costos.

Ahora, si miramos las varianzas estimadas para los efectos aleatorios, podemos ver que tanto PIB como Desempleo tienen varianzas considerables teniendo en cuenta la escala de los datos: 3.58 para Desempleo y 0.57 para PIB (la varianza del intercepto no parece ser muy notable). Esto nos demuestra lo que esperábamos: el componente aleatorio juega un papel muy importante en este caso ya que cada empresa tiene una dinámica muy particular y por ende el modelo le da un buen peso a los efectos aleatorios. Vemos entonces que el PIB y el desempleo afectan de manera muy diferente y particular a cada empresa. Además, pareciera verse que la forma en la que el desempleo afecta a cada empresa es mucho más notoria y particular y de mayor peso que el efecto que tiene el PIB en cada empresa (ya que la varianza del desempleo en los efectos aleatorios es mucho mayor a la del PIB).

Usamos el conjunto de entrenamiento para medir la capacidad de ajuste. En el conjunto de entrenamiento, el modelo de efectos mixtos propuesto se comporta así:

preds = predict(mod_me)

r2_model<-r2_score(preds, data_train_z$Costo.de.ventas_dif)
adj_r2_model<-adj_r2_score(preds, data_train_z$Costo.de.ventas_dif)
mse_model<-mean((data_train_z$Costo.de.ventas_dif - preds)^2)

#sub_tit = paste("R2", format(r2_model, digits=2, nsmall=2), 
#                "; R2_adj", format(adj_r2_model, digits=2, nsmall=2), 
#                sep = " ", collapse = NULL)

plot(lapply(preds, function(x) exp(x) - 1),
     lapply(data_train_z$Costo.de.ventas_dif, function(x) exp(x) - 1), 
     main='Modelo de Efectos Mixtos (variación porcentual anual)', #sub = sub_tit,
     xlab = 'Valor predicho', ylab = 'Valor Real')
abline(a=0, b=1, lwd = 2, col = 'red')

Medidas Resultados
r2 train 0.4016082
r2 adj train 0.3922583
mse train 0.1557839

Bajo esta aplicación de efectos mixtos es claro como se presenta una mejor distribución de los valores predichos contra los reales y éste es un comportamiento más natural y razonable en la industria en comparación al obtenido con el lineal general. Si bien hay unas compañías para las cuales su predicción en el conjunto de entrenamiento se aleja mucho del resto de las empresas, la predicción para la mayoría de las empresas permanece en escalas de entre -0.2 y 0.1 con una variedad muy interesante ya que se evidencia el efecto distribuido tanto del PIB como del desempleo en las variaciones de los costos de ventas en una serie de compañías que seguramente permanecen estables en el mercado. Respecto a las compañías que se hallan por fuera de este intervalo se tienen dos clases, por un lado aquellos cuyos valores reales son menores a cero son compañías cuyos costos de ventas disminuyeron más de lo que se esperaría, y por otro lado están los casos donde los valores reales fueron considerablemente superiores a los estimados, lo que refleja compañías cuyos costos de venta superaron por mucho las expectativas.

Se considera entonces que las predicciones obtenidas por el modelo en el conjunto de entrenamiento son razonables, ya que vemos que sus predicciones están todas entre -0.4 y 0.1, lo cual denota valores factibles en los cuales podría oscilar la variación porcentual de las empresas (tener reducciones máximas de 40% o incrementos de 10%), lo cuál es consistente con la media de -0.03 que habíamos calculado para la variable explicativa (en promedio, hay una disminución en cada año en los costos de ventas de las empresas tomadas). Además, se observa una gran variabilidad en las predicciones debido a la existencia de efectos aleatorios, algo que se esperaba ver en este modelo de efectos mixtos.

Adicionalmente, al ver las métricas obtenidas por este modelo en el conjunto de entrenamiento comprobamos que el ajuste realizado por este modelo es mucho más decente que el del modelo lineal general: tuvo un R2 de 0.4016 el cual representa un ajuste decente (el del lineal general era casi de 0). Además su R2 ajustado fue de a 0.3923 (decente). Por último, su MSE fue de 0.1558, lo cuál representa una mejora muy considerable en comparación al MSE del conjunto de entrenamiento del modelo lineal general que era de 0.2025. Estos resultados indican que, tal como se esperaba, el modelo de efectos mixtos logra realizar un ajuste mucho mejor al incluir los efectos particulares de cada empresa.

Por otro lado, veamos cómo se comporta el modelo de efectos mixtos en el conjunto de prueba. Recordemos que en este caso el modelo tendrá que hacer predicciones para individuos (empresas) que nunca había visto, por lo que deberá predecir usando solo el componente poblacional general, ya que no conoce el componente aleatorio particular que tendrían las empresas del conjunto de prueba.

preds = predict(mod_me, newdata = data_test_z,allow.new.levels=TRUE)

r2_model<-r2_score(preds, data_test_z$Costo.de.ventas_dif)
adj_r2_model<-adj_r2_score(preds, data_test_z$Costo.de.ventas_dif)
mse_model<-mean((data_test_z$Costo.de.ventas_dif - preds)^2)

sub_tit = paste("R2", format(r2_model, digits=2, nsmall=2), 
                "; R2_adj", format(adj_r2_model, digits=2, nsmall=2), 
                sep = " ", collapse = NULL)

plot(lapply(preds, function(x) exp(x) - 1),
     lapply(data_test_z$Costo.de.ventas_dif, function(x) exp(x) - 1),
     main='Modelo de Efectos Mixtos (variación porcentual anual)', xlab = 'Valor predicho', ylab = 'Valor Real')#, sub=sub_tit)
abline(a=0, b=1, lwd = 2, col = 'red')

Medidas Resultados
r2 test 0.0100776
r2 adj test -0.0218553
mse test 0.1805439

Desafortunadamente en el caso de prueba no se logra distinguir la tendencia natural en el comportamiento de la diferencia de los costos de venta y la escala se encuentra muy distante entre los valores predichos y los reales. En esta evaluación los valores predichos se ponderan para generar una tendencia casi constante, tal como era de esperar debido a que el modelo de efectos mixtos en este caso solo pudo utilizar el componente poblacional general. De esta forma, al usar solo el componente poblacional, el modelo usó solo los valores de las variables macro para predecir, y debido a que solo se tenían 3 valores únicos para cada variable macro, el modelo es capaz de generar solo 3 escenarios distintos de predicción, tal cómo se observa en la gráfica.

Adicionalmente, al ver las métricas obtenidas por este modelo en el conjunto de prueba corroboramos también que su capacidad predictiva es también pobre: el MSE es de 0.1805, mientras que el R2 y el R2 ajustado son muy bajos (0.0101 y -0.0219 respectivamente). Sin embargo, notemos que, a pesar de que el MSE dio muy levemente peor que el del modelo lineal general, el R2 y el R2 ajustado dieron mejores en este caso (y por una escala más considerable que la diferencia en los MSE). Esto daría indicios de que, a pesar de que el modelo de efectos mixtos aún no tiene una buena capacidad predictiva, fue capaz de separar un efecto poblacional y realizar un modelamiento un poco mejor que el modelo lineal general presentado anteriormente.

Observando ya las cifras predichas, vemos un comportamiento bastante similar al modelo lineal general en el conjunto de prueba. Esto era de esperar debido a que, cómo lo dijimos, el modelo de efectos mixtos al parecer depende mucho de los efectos aleatorios y en los efectos fijos los coeficientes eran similares al modelo lineal general (a pesar de que el modelo lineal general tenía solo el PIB y el de efectos mixtos tenía el PIB y el Desempleo, notar que las métricas en el conjunto de prueba fueron similares y el coeficiente del PIB y el intercepto de los efectos fijos es relativamente similar en ambos).

De acuerdo con lo dicho anteriormente, vemos que las predicciones en el conjunto de prueba oscilan desde -0.15 hasta -0.06, y como se dijo anteriormente para el modelo lineal general, es razonable que las predicciones oscilen entre estos valores. El modelo está prediciendo variaciones porcentuales anuales en los costos de las empresas de entre -6% y -15%, estos valores son bastante factiables ya que indican una disminución no muy grande de los costos, algo que es común que pase.

En resumen, luego de ver tanto el ajuste como la capacidad predictiva del modelo de efectos mixtos, vimos que el modelo de efectos mixtos logró obtener un ajuste mucho mejor que el lineal general. Sin embargo, al evaluar su capacidad predictiva en empresas que nunca antes había visto, como era de esperar sus resultados no fueron tan buenos y dependió completamente de su componente poblacional, el cuál presentaba variabilidad muy baja ya que solo se usaron 3 valores únicos en las variables explicativas. No obstante, el modelo logró detectar un leve componente poblacional y separarlo de los efectos aleatorios, y se aprecia una capacidad predictiva levemente superior a la del modelo lineal general (aunque sigue siendo muy baja).

Superficie poblacional del modelo de efectos mixtos

Miraremos la superficie poblacional del modelo, es decir, como cambia el valor predicho de una solución general desconocida de la población (para la cual solo se le aplicarían los efectos fijos, y no los aleatorios) al variar los valores de las dos variables explicativas que inciden sobre ella (desempleo y PIB). Se menciona que aquí se grafican las variables con la transformación logarítmica aplicada (no se hace la parte de devolver la salida a la escala de variación porcentual anual, sino que se dejan con el logaritmo).

Desempleo = seq(min(data_train_z$Desempleo),max(data_train_z$Desempleo), length.out = 50)
PIB = seq(min(data_train_z$PIB),max(data_train_z$PIB), length.out = 50)

data_surface = expand.grid(Desempleo,PIB)
names(data_surface) = c('Desempleo','PIB')
data_surface$NIT = 1000

z = predict(mod_me, newdata = data_surface,allow.new.levels=TRUE)
z = matrix(z, nrow = 50, ncol = 50)


custom_txt <- paste0("Desempleo: ", data_surface$Desempleo,
                    "\nPIB: ", data_surface$PIB, # correct break syntax
                    "\nCosto_de_ventas_dif: ", z) %>%
    matrix(50,50) # dim must match plotly's under-the-hood? matrix 


### Toca rotar la figura para poder ver el plano
fig <- plot_ly(x = PIB, y = Desempleo, z = z,    text = custom_txt,      hoverinfo = "text") %>% add_surface(colorbar=list(title='Poblacional'))%>%layout(scene = list(xaxis = list(title = 'PIB'), yaxis = list(title = 'Desempleo'),zaxis = list(title = 'Costo_de_ventas_dif')))
fig

Superficie del modelo para dos individuos conocidos conocidos por el modelo de efectos mixtos

Para visualizar mejor el funcionamiento del modelo de efectos mixtos, además de graficar la superficie poblacional general, graficaremos también la superficie de un par de invidiuos ya conocidos por el modelo. En este caso, como los individuos fueron usados por el modelo de efectos mixtos al entrenarse, se les aplica a dichos individuos un efecto aleatorio propio a cada uno de ellos (además del efecto fijo general que se aplica de la misma forma para todos). Por esta razón, veremos que las superficies de estos dos individuos son distintas entre sí y además son distintas de la superficie poblacional general.

Desempleo = seq(min(data_train_z$Desempleo),max(data_train_z$Desempleo), length.out = 50)
PIB = seq(min(data_train_z$PIB),max(data_train_z$PIB), length.out = 50)

data_surface = expand.grid(Desempleo,PIB)
names(data_surface) = c('Desempleo','PIB')
data_surface$NIT = empresas_train[1]

z2 = predict(mod_me, newdata = data_surface)
z2 = matrix(z2, nrow = 50, ncol = 50)

### Toca rotar la figura para poder ver el plano
data_surface$NIT = empresas_train[3]
z3 = predict(mod_me, newdata = data_surface)
z3 = matrix(z3, nrow = 50, ncol = 50)


fig %>% add_surface(z = ~z2, opacity = 0.98,colorscale = list(c(0,1),c("rgb(255,112,184)","rgb(120,0,64)")),colorbar=list(title=paste('NIT ', empresas_train[1]))) %>% add_surface(z = ~z3, opacity = 0.94,colorscale = list(c(0,1),c("rgb(0,0,0)","rgb(200,200,200)")),     colorbar=list( title=paste('NIT ', empresas_train[3]))) %>% layout(title="Poblacional (azul y verde) vs dos Individuos (uno rosa y uno negro)")

De las gráficas obtenidas al momento de mostrar la superficie del comportamiento de la variable de respuesta en base a la variación del PIB y el Desempleo, vemos que en particular las dos empresas que se muestran tienen un incremento porcentual anual en sus costos de venta mayor al poblacional general (ya que es claro que dichas superficies se ubican por encima de la superficie poblacional de color azul y verde, la cual nos indica el comportamiento poblacional general de la variable Costo_de_ventas_dif). Además, cabe resaltar que la empresa con NIT 800081030 (superficie negra) tiene una tendencia a tener un mayor aumento en la variable de respuesta en comparación a la empresa con NIT 800015615 (superficie rosa). Esto indicaría entonces que la empresa de construcción NIT 800081030 por lo general tiene mayor crecimiento (mayor aumento en sus costos) que la empresa con NIT 800015615, según lo captado por el modelo de efectos mixtos.

Conclusiones

Referencias

Estimación del esfuerzo (Anexo)

Según lo pedido en el enunciado del trabajo, entregamos aquí una estimación de las horas empleadas para desarrollar este trabajo. Las horas aquí mencionadas son horas totales de trabajo (es la suma de las horas empleadas por los 3 autores del reporte).

  1. Consolidación de la información : 14 horas

  2. Transformación de variables y análisis descriptivo: 18 horas

  3. Ajuste y validación de modelos: 18 horas

  4. Redacción del reporte: 23 horas